function [dist, bc] = gpdist(s1,s2)
% Compute distance between gp posteriors s1 and s2. Column one is the mean
% value, and Column two corresponds to the variance



% Split variables
up = s1(:,1); % mean s1
% up2 = s1(:,1).^2; % mean s1 squared
sp = s1(:,2); % var s1
sp2 = sp.^2; % var squared s1
uq = s2(:,1); % mean s2
% uq2 = s2(:,1).^2;
sq = s2(:,2);
sq2 = sq.^2; % var squared s2

% Compute Bhat. distance according to wikipedia definition:
%  - Non-vectorization:
% dist = zeros(size(s1,1),1);
% for i = 1 : length(dist)
%     dist(i) = .25*(log((sp(i)./sq(i))+(sq(i)./sp(i))+2) + ((up(i)-uq(i))^2)./(sp(i)+sq(i)));
% end
%  - Vectorisation
dist = .25*(log(.25*((sp./sq)+(sq./sp)+2)) + ((up-uq).^2)./(sp+sq));

bc = exp(-1*dist);
% Compute Bhat. COEFFICIENT according to Benjamin M. Marlin definition
%  - Non-vectorization:
% bc = zeros(size(s1,1),1);
% for i = 1 : length(bc)
%     ut = (up(i)/sp(i)) + (uq(i)/sq(i));
%     st = 1/((1/sp(i)) + (1/sq(i)));
%     ut2 = ut^2; st2 = st^2;
%     C = (sqrt(2)*st)/sqrt(sp(i)*sq(i));
%     
%     bc(i) = C*exp(-1*((up2(i)/sp(i))+(uq2(i)/sq(i))-(ut2/st2)));
% end
%  - Vectorisation
% ut = (up./sp) + (uq./sq);
% st = 1./((1./sp) + (1./sq));
% ut2 = ut.^2; st2 = st.^2;
% C = (sqrt(2)*st)./sqrt(sp.*sq);
% 
% bc2 = C.*exp(-1*((up2./sp2)+(uq2./sq2)-(ut2./st2)));



